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I.  INTRODUCTION 

With  the  recent  decline  in  the  threat  of  deep  water  submarine  forces,  growing  emphasis 
is  being  placed  on  operations  in  relatively  shallow  water.  In  contrast  to  open  ocean  models, 
shallow  water  models  of  sound  propagation  are  influenced  by  interaction  with  the  bottom. 
Variables  such  as  bottom  composition,  slope  angle,  sediment  density,  etc.,  make  acoustic 
predictions  difficult  at  best.  At  present  there  is  no  closed  analytic  solution  to  the  problem  of 
a  penetrable,  sloping  bottom  as  found  close  to  shorelines  or  along  the  continental  shelf. 
However,  computer  models  using  various  approximating  methods  are  in  use.  Computer 
models  are  limited  only  by  the  speed  and  memory  of  the  computer.  With  faster  computers 
being  developed,  more  computationally  intensive  models  become  practical.  These  models 
include  the  parabolic  equation  approximation,  image  theory,  and  adiabatic  normal  mode 
theory. 

A.    PARABOLIC  EQUATION 

By  replacing  the  Helmholz  equation  with  a  one-way  parabolic  equation  (PE),  an 
acoustic  field  can  be  generated  as  an  initial  value  problem.  The  PE  method  was  developed 
for  the  radar  community  and  first  applied  to  underwater  acoustics  by  Tappert  with 
restrictions  on  the  maximum  angular  aperture  at  the  source  [Ref  1].  Collins  finally 
eliminated  this  limitation  of  the  PE  model  with  a  higher  order  parabolic  equation  [Ref  2].  PE 
methods  by  Jensen  and  Kuperman  [Ref  3]  took  advantage  of  improved  computing  power 
and  showed  valid  results  when  compared  to  model  tank  experiments  [Ref  4]  for  upslope 
propagation.  This  model  showed  modal  cutoff  as  the  sound  traveled  upslope.  It  also  showed 
the  sound  propagating  into  the  bottom  at  the  modal  cutoff  depths.  Jensen  and  Tindle  used 
improved  numerical  methods  to  further  study  this  problem  [Ref  5]. 
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To  increase  accuracy,  Lee  and  McDaniel  broke  the  wedge  up  into  a  number  of  discrete 
sections  and  applied  the  normal  PE  methods  [Ref  6].  A  more  recent  method  of  rotated 
parabolic  equations  by  Collins  has  shown  good  results  [Ref  7].  Finally,  Fawcett  has 
developed  a  three-dimensional  computer  model  using  PE  methods  combined  with  Fast 
Fourier  Transform  (FFT)  algorithms  [Ref  8]. 

B.    METHOD  OF  IMAGES 

Similar  to  the  image  reflections  in  a  kaleidoscope,  this  method  combines  the 
contributions  of  the  virtual  sound  sources  to  calculate  the  pressure  field.  Macpherson  and 
Dainteth  proposed  a  phase  incoherent  model  for  upslope  propagation  in  1966  using  the 
method  of  images  with  some  success  [Ref  9].  Coppens  et.  al.  [Ref  10]  later  developed  a 
phase  coherent  model  to  calculate  the  pressure  amplitude  along  the  bottom  of  a  wedge  of 
water  overlying  a  fast  fluid  bottom  in  the  upslope  direction  using  the  incident  plane  wave 
reflection  coefficient  stated  by  Brekhovskikh  and  Lysenov  [Ref  1 1].  Baek  then  developed  a 
computer  model  of  pressure  throughout  the  water  over  the  fast  bottom  [Ref  12].  However, 
this  model  was  limited  to  a  fast  bottom  due  to  the  critical  grazing  angle  of  the  sound  to  the 
bottom.  LeSesne  used  the  same  model  as  a  basis  for  a  three  dimensional  model,  which 
produced  pressure  fields  in  the  cross  slope  case  [Ref  13].  He  further  validated  his  computer 
results  with  experimental  results  from  a  model  tank.  Kaswandi  then  developed  the  model  for 
a  slow  bottom  in  the  dovvnslope  direction  [Ref  14].  Concurrently,  Li  introduced  a  model 
which  applies  to  both  fast  and  slow  bottoms  [Ref  15].  Nassopoulos  then  took  all  the 
previous  models  under  development  at  the  Naval  Postgraduate  School  and  combined  them 
in  the  computer  program  URTEXT  [Ref  16]  Independent  of  this  research,  Buckingham  and 
Deane  developed  a  similar  model,  which  produced  a  transmission  loss  plot  as  a  function  of 
depth  and  range  from  the  source  [Ref  17] . 


The  analytic  development  of  the  method  is  an  ongoing  process.  Nassopoulos 
implemented  the  use  of  image  doublets  [Ref  16].  Livingood  then  developed  the  image 
doublets  to  the  three  dimensional  case  [Ref  18].  Presently,  work  on  Taylor  expansions  of 
sets  of  four  images  is  being  done  by  Joyce  [Ref  19]. 

Tests  done  by  Kim  [Ref  20]  and  by  LeSesne  are  cited  to  prove  the  validity  of  the  image 
model.  LeSesne  compared  his  model  with  experimental  conditions,  while  Kim  compared  his 
to  a  parabolic  equation  model  developed  by  Jaeger  [Ref  21]  which  used  an  implicit  finite- 
difference  solution.  This  turned  out  to  be  a  favorable  comparison. 

While  validated  in  comparisons  with  both  other  methods  and  with  experiment,  the 
method  of  images  still  has  some  flaws.  It  is  very  computer  intensive,  especially  as  the  angle 
of  the  bottom  becomes  very  shallow.  Improvements  include  truncation  of  images  that  do  not 
significantly  contribute  to  the  total  pressure  field. 

C.    NORMAL  MODES 

Normal  mode  theory  is  a  range  independent  approach.  However,  for  the  wedge 
problem  a  range  dependent  model  is  needed.  Pierce  showed  that  an  approximation  of  the 
normal  modes  can  be  used  by  performing  an  adiabatic  separation  of  depth  and  range  in  the 
wave  equation  [Ref  22].  Graves,  Nagl,  Uberall,  and  Zarur  then  applied  this  method  to  the 
wedge  problem  using  isospeed  water  over  a  rigid  bottom  [Ref  23].  This  method  achieved 
good  agreement  at  large  angles  when  compared  to  experimental  values  [Ref  24].  This 
approach  was  limited  to  large  angles  and  a  perfectly  reflecting  bottom.  Buckingham 
developed  a  solution  of  the  penetrable  bottom.  In  this  method  he  first  calculates  an  effective 
pressure  release  bottom  below  the  actual  penetrable  bottom.  The  depth  of  this  bottom  is  a 
function  of  the  density  and  speed  of  sound  of  the  actual  bottom.  He  then  applied  the  sum  of 
the  normal  modes  to  this  new  bottom  [Ref  25]. 


A  problem  of  the  adiabatic  normal  mode  approximation  is  that  it  does  not  adequately 
explain  the  transition  to  evanescent  modes  at  cutoff  depths  [Ref  3].  By  combining  PE  with 
adiabatic  mode  theory',  a  critical  depth  function  was  developed  which  agreed  with  the  results 
of  Jensen  and  Kuperman  [Ref  26]. 

Another  method  is  to  combine  ray  acoustics  with  local  modes  [Ref  27],  which  leads  to 
the  development  of  intrinsic  wave  modes  [Ref  28].  Also  using  beam  displacement  with  ray 
theory  has  shown  agreement  with  the  two-way  coupled  mode  solution  in  two  dimensions 
over  a  penetrable  lossy  bottom  [Ref  29]. 

To  take  into  account  coupling  between  the  modes,  a  version  of  stepwise  depth 
variations  was  applied  to  normal  mode  theory  [Ref  30].  This  method  was  used  first  by  Chiu, 
Desaubies,  and  Miller  [Ref  31],  then  by  Chiu  and  Ehret  into  the  development  of  a  model  of 
coupled  normal  modes  [Ref  32].  The  three  dimensional  analysis  by  Sagos  [Ref  33]  is  the 
extent  of  published  information  on  coupled  normal  modes,  but  improvements  to  the  theory 
and  the  model  are  still  being  pursued. 


II.  DEVELOPMENT 

A.    IMAGE  THEORY 

In  1978,  Coppens,  Sanders,  Iounnou,  and  Kawamura  predicted  and  measured  the 
pressure  field  and  phase  up  slope  along  a  wedge  shaped  bottom  [Ref  10].  LeSesne 
implemented  a  model  that  was  not  dependent  on  up  or  down  slope  direction  [Ref  13]. 
Further  research  cited  earlier  developed  cross  slope  programs,  giving  three  dimensional 
capability  to  the  model.  The  theory  behind  these  models  is  the  method  of  images. 

Coppens  et.  al.  provides  the  complete  explanation  of  the  theory  but  a  summary  is  in 
order.  In  transmission  within  a  waveguide  there  is  a  number  of  interactions  of  acoustic  rays 
with  the  guide  boundaries  [Ref  34,  p.  427].  In  image  theory  each  ray  path  is  replaced  by  an 
image  of  the  source,  the  distance  of  the  source  from  the  receiver  being  the  total  distance  of 
the  actual  acoustic  ray.  For  a  sloping  bottom  the  image  is  placed  a  distance  from  the 
reflecting  boundary  equal  to  the  distance  of  the  source  from  the  boundary  perpendicular  to 
that  boundary.  Figure  1  illustrates  this  geometry. 

The  phase  of  the  image  is  180°  out  of  phase  from  the  source  if  there  is  an  interaction 
with  a  pressure  release  surface,  as  in  the  case  of  the  air-water  interface,  and  in  phase  for  a 
rigid  surface— a  hard  bottom,  for  example.  Each  interaction  of  the  sound  with  the 
boundaries  implies  the  existence  of  a  new  image.  The  final  pressure  at  the  receiver  can  be 
calculated  by  taking  the  sum  of  the  contributions  from  each  image. 


Figure  1  -  Two  Dimensional  Wedge  Geometry 

For  a  given  bottom  slope  angle  p,  the  number  of  images  N  in  each  of  the  upper  and 
lower  half  spaces  is  given  by 

where  (3  is  an  integer  factor  of  180°.  Each  image  is  numbered  from  the  source,  with  the 
source  being  image  number  1  for  the  upper  half  space.  For  the  lower  half  space  the  first 
image  below  the  bottom  is  image  number  1. 

Next,  the  grazing  angle  6„  of  the  //th  image  above  or  below  the  bottom  is  calculated 
from 


8n  =  «p-Y  for  n  odd  (2a) 

8n  =  (n  -  l)p  +  y  for  n  even  (2b) 

where  y  is  the  angle  of  the  source  from  the  surface. 

The  range  Rn  from  each  image  to  the  receiver  is  then  calculated  using 

^=Vri2+r2-2r,r2cos(8n-p  +  6) 


for  the  upper  series  of  images  and 

K  =^f+'-2-2r1r2cos(e„+p-6) 


(3a) 


(3b) 


for  the  lower  series.  In  these  equations,  r\  is  the  range  from  the  wedge  apex  to  the  source,  ri 
is  the  range  from  the  apex  to  the  receiver,  and  6  is  the  angle  from  the  surface  to  the  receiver. 

For  each  interaction  with  the  surface,  the  reflection  coefficient  is  -1.  Each  interaction 
with  the  bottom  requires  calculation  of  the  reflection  coefficient.  This  coefficient  is  a 
function  of  the  speed  of  sound  in  each  medium,  the  density  of  each  medium,  and  the  grazing 
angle  of  the  ray  equivalent  on  the  bottom.  It  is  sufficient  to  find  the  sine  of  this  angle  6,„„. 
For  this 

rx sin( 6W  -2wp)  +  r2sin[(2w-l)p+6] 


sin6„m  = 

nm 


R 


n 


(4a) 


for  the  upper  images  and 


sine        rlSin(en-2mP)+r2sin[(2;n  +  l)p-6] 
for  the  lower  images  [Ref  16].  Figure  2  illustrates. 
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Figure  2  -  Incident  Angle  6„m  Calculation  (Upper  Half  Space) 

The  reflection  coefficient  for  each  interaction  R^,  then  is  given  by 


H-nm  ~ 


^smQnm  +  ^[^r^+a)y2-j^[^77-af 


(5) 


where  pi,  02  are  densities  of  the  respective  media,  j  =V-1,  and  a  and  b  are  given  by 

2 


<7  = 


\ 


<V 


-cos  0. 


b  =  2 


cx  i     a 


\  c2/    A^ 


(6a) 


(6b) 


For  the  above,  ch  C2  are  the  speeds  of  sound  in  the  respective  media  and  alki  is  a 
measure  of  the  bottom  loss  coefficient.  [Ref  1 1].  We  further  define  kn  as  the  wave  number 
in  the  wth  medium,  where 

K  -  °/c  (6c) 
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The  pressure  for  each  half  space  is  given  by  summing  the  contribution  for  the  upper 
and  lower  images  respectively 


n-l  **„  m-1 

where  I  Tr  =  1  when  n  =  1,2;  in  (7a)  and 

P^S^expC-ytAX-l^rK 


(7a) 


■l"fl 


(7b) 


The  number  of  interactions  m  with  the  bottom  is  related  to  the  image  number,  so  the  limit  M 

- 1  for  the  Arth  image.  For  convenience  we  have  omitted  the  factor 


is  the  integer  of 

eim'  which  would  otherwise  appear  in  (7a)  and  (7b). 

Two  special  cases  exist  in  the  upper  half  space.  For  the  case  of  the  direct  path  of  the 
source  to  the  receiver  (n=\),  R  is  1.  In  the  case  of  a  single  reflection  off  the  surface  and  no 
reflections  with  the  bottom  (n=2),  R  is  again  1.  The  phase  inversion  of  the  ray  is  computed 
by  the  (-1)^  term  of  (7a). 

Since  the  model  is  frequency  independent,  we  will  scale  the  ranges  to  a  distance  Xc 
from  the  apex  where  the  bottom  depth  is  equal  to  the  lowest  normal  mode  cutoff  depth  when 
the  bottom  is  fast  (i.e.,  when  the  sound  speed  ratio  is  less  than  1).  This  cutoff  distance  is 
calculated  by 


*£.- 


it 


2sin0ctanj3 


6.  =  cos 


-i 


I '    \ 


\C2/ 


(8a) 


For  a  slow  bottom  an  analogous  scaling,  convenient  for  the  computer  program  is 


k{Xc  = 


ji 


0.  =cos 


\C1/ 


(8b) 


2tan05tan/S 
For  the  total  pressure  at  the  receiver,  the  upper  and  lower  contributions  are  combined, 


thus 


p  -  Pu  +  p/  (9) 

B.     COUPLED  MODE  THEORY 

Normal  mode  theory  describes  sound  propagation  as  a  collection  of  eigenfunctions, 
called  normal  modes,  which  are  determined  by  the  depth  of  the  sound  channel  boundaries 
and  the  source  depth.  Coupling  between  the  modes  (energy  transfer)  in  the  horizontal 
direction  occurs  due  to  changes  in  sound  speed,  density,  or  bathymetry.  The  Chiu-Ehret 
model  [Ref  32]  determines  the  pressure  field  using  numerical  evaluations.  Sagos  [Ref  33] 
modified  this  model  for  three  dimensions. 

Chiu  [Ref  35]  then  developed  a  model  in  MATLAB™  called  BBCM.  In  this  model  the 
spectral  pressure  at  the  receiver  is  decomposed  by  an  FFT  algorithm.  At  each  point  in  the 
pressure  field  and  at  each  frequency,  the  amplitude  envelope  for  each  existing  normal  mode 
is  calculated,  then  coupled  to  determine  the  total  pressure  by 


^z,/)  =  /5(/)|^^(r;/)Zn(z;r,/)exp|/(j^2^n(r1;/Vr-2^) 


dt      (10) 


where  S(f)  is  the  source  amplitude  spectrum,  Z„  is  a  linear  combination  of  the  normal 
modes,  and  U^r;f)  is  the  amplitude  slowly  varying  modulation  complex  envelope.  It  is  in 
this  calculation  that  the  coupling  of  the  normal  modes  occurs. 

In  the  MATLAB™  version,  it  has  several  advantages: 

•  Can  handle  both  discrete  and  broad  band  signals  defined  by  S(f), 

•  The  bottom  can  be  varied— not  limited  to  a  smooth  slope, 

•  FFT  algorithms  are  extremely  fast.  Pressure  fields  are  calculated  in  a  fraction  of 
the  time  that  would  be  taken  by  URTEXT. 

BBCM  uses  a  rigid  bottom  for  the  waveguide,  while  in  the  upcoming  comparison  no 

such  bottom  is  assumed.  Therefore,  the  effect  of  this  rigid  bottom  on  the  sound  in  the  water 

column  must  be  examined. 
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C.    THE  TWO  LAYER  WAVE  GUIDE 

Figure  3  shows  a  sound  velocity  profile  of  a  hypothetical  wave  guide.  In  this  guide  the 
boundary  conditions  are  a  pressure  release  surface  and  a  rigid  bottom,  with  a  discontinuous 
change  in  sound  speed  at  depth  L 


z  =  0 

CO 

CO 
Q 

z=//  , 

Speed  of  sound 

C=C, 

c  =  c2 

Figure  3  -  Sound  Velocity  Profile  for  Comparison 

Normal  mode  theory  produces  an  exact  solution  to  the  normalized  Helmholtz  equation 
for  sound  pressure  P.  This  equation,  given  in  cylindrical  coordinates  for  a  point  source,  is 


r  dr 


dr 


p  +  4p  + 

dz1 


O) 


c(z) 


P  =  -b(r)b{z-z0)e 
r 


JWl 


(11) 


and  has  a  solution  of 

P(r,z,t)  =  -j^'^ziz0)Zn{z)[l%\Knr)]  (12) 

n 

where  Zn(z)  is  the  depth  dependent  normal  mode,  z„  is  the  depth  of  the  source,  r  is  the 
range  from  the  source,  Kn  is  the  horizontal  propagation  constant  of  the  nth  normal  mode, 
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and  f^2)is  the  Hankel  function  where  #J2)  -JQ  -  jY0  [Ref  34,  p.  431].  The  function  Zn(z) 

is  a  solution  to  the  homogeneous  Helmholtz  equation 

d2 
-^rZn+k*Z=0  (13) 

dz  z 

First  we  examine  the  case  where  H-  oo.  Zn(z)  must  satisfy  the  boundary  conditions  of 

Z(0)=0,  continuity  of  particle  velocity  and  pressure  across  the  interface  at  z  -  L,  and  Zn(z) 

must  approach  zero  as  z  goes  to  infinity. 

For  the  above  conditions,  we  guess  at  two  functions  Z,„  and  Z2n  to  represent  Zn(z), 

Zln(z)'Ansm(klvlz)  0<z<L  (14) 

and 

Z2n(z)  =  Bne-p-{z-L)  L<z  (15) 

where  kUn  is  the  vertical  wave  number  of  the  wth  depth  dependent  normal  mode  in  the  upper 
layer  and  p„  is  the  vertical  propagation  number  of  the  exponential  function  in  the  lower 
layer.  The  numbers  klzn  and  p„  are  related  by  the  equation 


2 
K_  = 


^ 


-*L- 


/  v 

*     +ti  (16) 

\C2/ 


\Cl/ 

Evaluation  of  (14)  at  the  surface  proves  that  it  is  a  valid  solution,  but  the  behavior  of  the 
functions  at  the  interface  of  z  =  L  needs  to  be  determined.  At  the  interface,  continuity  of 
sound  pressure  and  normal  particle  velocity  must  hold.  Therefore, 

AM^nL)=Bn  (17) 

and 

4fc„cos(*laiL)--  Bfiu  (18) 

Combining  (17)  and  (18)  and  manipulating  gives  the  transcendental  equation 
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PnL 


(19) 


which,  by  defining 


y  =  KzL 


(20) 


and 


2         r2 

a    =  L 


w 


\ci/ 


f  ,.A 


CO 


\C2/ 


(21) 


yields 


tony  =- 


v^ 


r 


(22) 


This  equation  can  be  used  to  solve  for  kUn,  and  this  yields  p„and  k„,  when  y  <  a. 
Normalization  of  Zn(z)  then  produces  the  terms  An  and  B„,  giving  the  final  solution  to  ( 13). 

This  solution  corresponds  to  an  oscillatory  function  above  L  overlying  a  function  which 
decays  exponentially  as  z  increases.  An  example,  Zj(z),  at  a  frequency  somewhat  above 
cutoff,  is  shown  in  Figure  4.  As  can  be  seen,  the  sound  energy  is  trapped  in  the  upper  layer 
with  some  energy  leaking  into  the  bottom. 
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Figure  4  -  Fluid  Layer  Overlying  an  Infinite  Bottom 
The  cutoff  frequency  of  each  normal  mode  can  be  determined  at  each  point  where  an 
approaches  an  asymptote  of  tan  y  in  (22),  or  where 


ji  3ji        (2w  -  1)ji 
a"  =  2'~,"',        2~ 


«  =1,2. 3. 


(23) 


Substitution  of  (21)  into  (23)  and  manipulation  yield  a  cutoff  angular  frequency  of 
normal  mode  n, 

(2w  - 1)^',^ 


(o.  = 


2L 


I^i) 


(24) 
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We  now  examine  the  case  where  there  is  a  rigid  bottom  at  z  -  H.  The  solutions  of  the 

dZ(z) 

Helmholtz  equation  now  must  fulfill  the  additional  boundary  condition  =  0 

dz    k.H 

For  the  revised  boundary  conditions,  we  choose 

Z,n{z)  =  A, sin(*lwz)  0<z<L  (25) 

and 

ZJn(z)  =  fi„cosh[pn(//-z)]  L<z<H  (26) 

The  hyperbolic  cosine  was  chosen  in  analogy  with  the  exponential  choice  in  the  previous 
case.  It  satisfies  the  boundary  conditions  at  z  =  H. 

Next  is  evaluation  of  (25)  and  (26)  at  the  interface  of  z  =  L.  At  the  interface, 

^sin(*u„L)  =  Bn  cosh[pn(//  -  L)]  (27) 

and 

AJclznoos(kiznL)  =  -fl„p„  sinh[pn(//  -  Lj\  (28) 

Combining  (27)  and  (28)  gives 

which,  by  substitution  ofy  and  a  yields 


tany 


ycoth[(%-l)V7T7' 


4^? 


(30) 


This  equation  can  be  used  to  solve  for  kUn  and  then  (3„and  Kn  can  be  obtained.  When  y<a 
the  solution  is  similar  to  that  for  the  infinite  bottom  condition  in  that  it  corresponds  to  an 
oscillatory  function  above  L  and  decays  as  z  increases  from  L  to  H.  The  sound  energy  in 
this  case  is  trapped  in  the  upper  layer  with  some  leakage  into  the  lower  layer  (Figure  5). 
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Note  that,  from  (30),  as  //decreases  toL,  y  decreases  to  (n  -  \J2)ji  for/z=  1, 2, 3,... 
and  the  normal  modes  will  tend  to  approach  those  for  a  singe  layer  waveguide  with  a  rigid 
bottom  and  a  pressure  release  surface. 


0 
-20 
-40 
-60 
-80 


£     -100 

Q 

-120 


-140 
-160 
-180 
-200 


1 

FS^"      1 

>. 

J 

-\i 

L=150m 

i 

\        jH=200  ra 

1  -0.5  0  0.5  1 

Normalized  Amplitude  Z 


Figure  5  -  Two  Layers  Overlying  a  Rigid  Sub-bottom 

The  propagation  of  the  normal  mode  exists  in  three  regimes  dependent  upon  klzn.  We 
determine  p„  in  terms  of  klzn  from  ( 16)  to  be 


fi- 

7    \ 

U) 

1 

[•] 

2" 

_KCU 

\C2) 

— -k 


Izn 


(31) 


16 


(m\ 


If*.\.< 


(I) 


\cl/ 


'at 


\C2/ 


,  this  is  equivalent  to  y  <  a,  and  therefore  (25)  through  (30)  hold.  If 


y  >  a  then  |3W  becomes  imaginary  and  Z^z)  becomes 

ZlH(z)—Bnca{kJLH-z)] 

where  kUn  is  defined  by 


(32) 


*2    =  fc2   - 


fj\ 


CO 


\C1/ 


f  ,A 


CO 


\C2/ 


(33) 


and  is  related  to  Kn  by 


*-\3  -'» 


(34) 


By  applying  the  boundary  contains  to  (25)  and  (32)  and  making  the  appropriate 
substitutions  for  y  and  a  we  produce 


tony  =- 


ycot 


Wl-W1* 


Vyw 


(35) 


This  equation  describes  a  normal  mode  trapped  in  the  waveguide  with  a  pressure  release 
surface  and  a  rigid  bottom.  How  that  mode  propagates  is  dependent  upon  Kn . 


/,.A2        /        "2 


From  (16),  if 


CO 


U, 


CO 


\H 


CO 


—     <  kizn  <    —    '  men  K„  is  real>  ar>d  the  mode  propagates 


■\ 


horizontally  in  the  overall  waveguide.  If    —     <  £2     tcn  becomes  imaginary.  The  normal 

mode  thus  decays  as  it  propagates  horizontally.  Figure  6  shows  the  vertical  component  of 
the  normal  mode  which  satisfies  (35). 
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Figure  6  -  Two  Layers  Overlying  a  Rigid  Sub-bottom 

The  condition  where  p„  becomes  imaginary  is  the  major  difference  between  the  infinite 
bottom  and  the  rigid,  finite  bottom.  In  the  infinite  bottom  case,  decreasing  the  layer  depth  L 
decreases  the  number  of  modes  which  can  be  propagated.  However,  the  finite  bottom  adds  a 
propagation  mode  throughout  waveguide  which  must  be  corrected  if  an  accurate  comparison 
is  to  be  made. 

Another  difficulty  encountered  is  that  the  above  discussion  is  valid  only  when  all 
interfaces  are  parallel.  With  parallel  layers,  the  grazing  angle  6  is  easily  computed.  With  a 
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sloping  interface  L=f(r)  this  angle  changes  with  range.  It  can  be  approximated  numerically, 
but  no  analytic  solution  exists. 

To  extend  to  a  sloping  interface  in  a  waveguide  consisting  of  two  layers  with  a  rigid 
bottom  of  depth  //and  fluid  interface  at  depth  L,  all  the  normal  modes  with  cutoff 
frequencies  below  the  frequency  of  the  source  are  excited  and  propagate  according  to  (25) 
and  (26).  As  the  sound  goes  upslope,  each  normal  mode  approaches  its  first  transition  depth 
where  jSn  =  0  and  the  pressure  amplitude  of  the  exponential  function  becomes  uniform  over 

the  entire  lower  layer  (Figure  7).  After  the  transition,  the  normal  mode  now  behaves 
sinusoidally  with  depth  in  both  layers  (Figure  6).  As  L  continues  to  decrease,  the  mode 
approaches  the  second  transition  point,  where  Kn  becomes  imaginary  and  the  mode 

becomes  evanescent  in  the  horizontal  plane. 
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Figure  7  -  Fluid  Layer  over  a  Rigid  Bottom  at  First  Transition 
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III.  PROCEDURES  AND  RESULTS 

A.    PROCEDURES 

For  this  research,  the  image  model  (URTEXT)  and  the  coupled  mode  model  (BBCM) 
were  compared  for  the  upslope  case.  Both  models  were  run  using  MATLAB™,  with 
URTEXT  being  translated  from  FORTRAN.  Validation  of  this  translation  appears  in 
Appendix  A. 

The  parameters  chosen  were  the  benchmark  wedge  problem  of  the  Journal  of  the 
Acoustic  Society  of  America  [Ref  36]  (Figure  8).  The  speed  of  sound  in  the  water  is 
uniform  at  1500  m/s  and  the  standard  seawater  density  of  1.021  g/cm3  is  used.  A  point 
source  transmitting  a  continuous  wave  at  25  Hz  is  located  4  km  from  the  apex  of  the  wedge 
and  100  m  below  the  surface.  The  bottom  slopes  at  2.86°.  The  composition  of  the  bottom 
will  be  varied  to  test  different  scenarios. 
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Figure  8  -  JASA  Benchmark  Wedge 
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Three  cases  were  run  for  comparisons.  Scenario  1  uses  a  pressure  release  bottom, 
Scenario  2  has  a  fast  bottom  with  no  absorption  loss,  and  Scenario  3  uses  the  same  bottom 
with  some  absorption.  Parameters  for  the  environments  are  given  in  Table  1. 
TABLE  1  -  ENVIRONMENTAL  PARAMETERS 


Scenario  1 

Scenario  2 

Scenario  3 

cj  (m/s) 

1500 

1500 

1500 

C2  (m/s) 

343 

1700 

1700 

Pi(g/cm3) 

1.021 

1.021 

1.021 

p2(g/cm3) 

0.00121 

1.500 

1.500 

a  (dB/X) 

0.00 

0.00 

0.5 

One  difficulty  was  the  difference  in  measured  values  between  the  two  models.  BBCM 
uses  SI  units  and  produces  transmission  loss  (among  other  values)  in  dB.  URTEXT  uses 
dimensionless  parameters  scaled  to  the  dump  distance  determined  by  (8a)  and  produces  a 
dimensionless  pressure  amplitude.  This  can  be  converted  to  dB  by 


TL 


-20,og(l% 


(36) 


where  TL  is  the  transmission  loss  in  dB,  P  is  the  dimensionless  complex  pressure,  and  Xc  is 
the  dimensionless  dump  distance.  For  each  scenario,  it  was  necessary  to  convert  the  SI  units 
into  the  dimensionless  input  units  for  URTEXT.  Actual  numbers  will  be  given  later. 

BBCM  uses  input  files  which  model  the  environment.  These  files  define  the  overall 
waveguide  (water  column  and  underlying  bottom)  by  describing  the  sound  velocity  profile, 
density  profile,  receiver  ranges,  receiver  depths,  attenuation  of  each  normal  mode,  and  source 
amplitude  spectrum  [Ref  35].  Depending  upon  the  degree  of  resolution  desired,  these  input 
files  can  get  quite  large.  A  program,  SCENEGEN,  was  written  to  develop  these  input  files 
specifically  for  the  two  dimensional  wedge. 
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Another  problem  of  BBCM  concern  the  rigid  bottom  of  the  overall  waveguide.  As  seen 
from  the  development  in  Chapter  II  this  rigid  waveguide  floor  affects  the  sound  in  the  water. 
Depending  on  the  overall  depth  of  the  waveguide,  free  propagating  normal  modes  of  the 
entire  waveguide  are  excited  which  are  coupled  with  the  modes  which  are  trapped  only  in  the 
water  layer.  These  free  normal  modes  must  be  identified  and  eliminated  for  accurate  results. 

This  could  be  done  by  applying  a  range  dependent  attenuation  coefficient  for  each 
normal  mode,  so  the  mode  is  attenuated  by  the  bottom  after  it  reaches  its  first  transition 
point.  However,  BBCM  does  not  support  this.  The  attenuation  vector  that  is  entered  in  the 
model  is  range  independent.  Therefore,  the  effects  of  the  rigid  sub-bottom  cannot  be 
eliminated  for  this  problem. 

B.    SCENARIO  1 -PRESSURE  RELEASE  BOTTOM 

A  pressure  release  bottom  was  simulated  for  both  models.  Since  both  models  are 
incapable  of  running  if  either  C2  or  p2=  0,  an  air- water  interface  was  deemed  sufficient  for 
the  bottom.  Inputs  for  URTEXT  are  given  in  Table  2.  As  the  sound  goes  closer  to  the  apex, 
adiabatic  theory  predicts  less  pressure  as  cutoff  depth  for  each  successive  normal  mode  is 
reached.  The  actual  outputs  are  seen  in  Figure  9. 

TABLE  2  -  URTEXT  INPUTS  FOR  COMPARISONS 


Scenario  1 

Scenario  2 

Scenario  3 

p 

3.0° 

2.86° 

2.86° 

Y 

1.43° 

1.43° 

1.43° 

y 

4.412 

0.8824 

0.8824 

y 

/P2 

850 

0.6667 

0.6667 

Xc 

69.87  m 

638  m 

638  m 

r\ 

57.33  Xc 

6.27^ 

6.21  Xc 

%2 

0.0001  dB 

0.0001  dB 

0.07963  dB 
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For  the  inputs  from  Table  2,  normal  mode  theory  for  a  rigid  bottom  predicts  seven 
modes  will  be  excited  by  the  source.  But  since  the  source  is  located  at  the  pressure  node  of 
the  even  numbered  modes,  these  modes  are  not  excited.  Therefor,  only  four  modes  are 
excited  and  cutoff  points  for  these  modes  should  be  observed.  The  URTEXT  results  show 
these  predicted  characteristics,  with  cutoff  points  for  the  first  mode  at  3400  m,  the  third 
mode  at  2200  m,  and  the  fifth  mode  at  1000  m.  Cutoff  for  the  seventh  mode  cannot  be  seen 
due  to  the  high  density  of  contour  lines  in  Figure  9.  Further  information  is  seen  in  the 
transmission  loss  plots  of  Figures  10  through  13. 

BBCM  did  not  produce  usable  data  for  this  scenario.  This  is  not  surprising,  when  one 
considers  the  original  purpose  of  the  model.  BBCM  was  written  for  acoustic  tomography 
calculations.  As  such,  the  pressure  release  case  is  not  considered,  since  there  are  no  actual 
pressure  release  bottoms  in  the  ocean.  The  model  ends  up  taking  the  sound  from  the  water 
column  and  trapping  it  in  the  underlying  bottom. 
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Figure  9  -  Scenario  1  -  URTEXT  Transmission  Loss  in  dB  re  1  m 
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Figure  10  -  Scenario  1  -  30  m  Receiver 
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Figure  11  -  Scenario  1  -  50  m  Receiver  Depth 
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Figure  12  -  Scenario  1  -  100  m  Receiver  Depth 
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Figure  13  -  Scenario  1  -  150  m  Receiver  Depth 


26 


C.    SCENARIO  TWO -FAST,  LOSSLESS  BOTTOM 

In  this  scenario  a  more  realistic  environment  was  used.  Inputs  were  as  seen  in  Table  1. 
This  produced  the  inputs  for  URTEXT  as  given  in  Table  2.  Figures  14  through  17  are 
transmission  loss  profiles  of  both  URTEXT  and  BBCM  within  the  water  column. 

Only  the  water  column  is  considered  as  URTEXT  does  not  calculate  pressures  in  the 
underlying  bottom.  However,  areas  of  significant  difference  can  be  seen.  These  are 
primarily  due  to  the  effect  of  the  rigid  sub-bottom.  Areas  where  BBCM  show  an  increase  in 
transmission  loss  are  likely  due  to  a  normal  mode  going  through  transition,  while  areas 
where  there  is  significant  decrease  in  transmission  loss  are  possibly  due  to  modal 
interference. 
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Figure  14  -  Scenario  2  (Lossless  Bottom)  -  30  m  Receiver  Depth 
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Figure  15  -  Scenario  2  (Lossless  Bottom)  -  50  m  Receiver  Depth 
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Figure  16  -  Scenario  2  (Lossless  Bottom)  -  100  m  Receiver  Depth 
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Figure  17  -  Scenario  2  (Lossless  Bottom)  -  150  m  Receiver  Depth 

D.    SCENARIO  3  -  LOSSY  BOTTOM 

The  environment  for  this  run  was  the  same  as  for  Scenario  2,  except  a  loss  term  of  0.5 
dB/X  was  added  for  the  bottom.  This  converted  to  0.07963  dB  for  the  loss  term  in 
URTEXT.  Figures  18  through  21  show  the  results  of  this  run. 

For  this  scenario,  another  difficulty  was  encountered.  While  the  benchmark  calls  for 
attenuation  by  the  bottom  for  the  source  frequency,  that  attenuation  must  be  applied  to  the 
normal  modes  of  interest.  Therefore,  a  range  dependent  attenuation  coefficient  must  be 
determined.  Since  no  such  coefficient  can  be  applied,  the  results  developed  are  rough 
approximations  at  best  near  or  around  modal  cutoff. 

BBCM  runs  were  not  accurate  but  are  estimates.  BBCM  was  run  with  three  normal 
modes.  Application  of  the  attenuation  term  is  a  difficult  exercise,  first  because  it  needs  to  be 
range  dependent  to  more  accurately  eliminate  the  rigid  sub-bottom,  and  second  because  the 
term  must  be  distributed  among  the  modes.  Further  research  in  this  problem  is  necessary. 


29 


Again,  good  correlation  is  noted  between  the  models.  The  addition  of  the  attenuation 
term  has  no  great  effect  on  the  sound  in  the  water  column.  However,  there  is  a  significant 
drop-off  at  the  end  of  the  BBCM  graphs.  These  are  the  points  where  the  receivers  are  now 
in  the  bottom.  Here  the  attenuation  is  significant. 

Another  significant  feature  is  in  Figure  18  at  150  m  range.  Here  a  dip  is  seen  in 
URTEXT  which  is  not  seen  in  BBCM.  At  this  relatively  short  range  surface  interference  is 
noted.  URTEXT  picks  this  up  (see  Appendix  A)  but  BBCM  does  not.  Again,  the  original 
purpose  of  BBCM  is  to  make  predictions  at  relatively  long  ranges,  significantly  greater  than 
the  influence  of  surface  interference  [Ref  35].  At  these  longer  ranges,  BBCM  is  a  robust 
model. 
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Figure  18  -  Scenario  3  (Lossy  Bottom)  -  30  m  Receiver  Depth 
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Figure  21  -  Scenario  3  (Lossy  Bottom)  -  150  m  Receiver  Depth 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 

Figures  22  through  25  show  the  results  of  Jensen  and  Ferla  using  a  PE  method 
compared  to  the  results  of  URTEXT  [Ref  36]. 
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Figure  22  -  Scenario  2  -  30  m  comparison  of  PE  and  URTEXT 
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Figure  23  -  Scenario  2  -  150  m  Comparison  of  PE  and  URTEXT 
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Figure  24  -  Scenario  3  -  30  m  Comparison  of  PE  and  URTEXT 


34 


-30 


70  —   I  I  I  I      I  I  I  I     I  I  I  I      I  I  I  I   TTT 


"\ 


URTEXT 
PE 


0      100     200     300     400    500     600     700     800     900    1000 
Range  from  Source  (m) 


Figure  25  -  Scenario  3  -  150  m  Comparison  of  PE  and  URTEXT 

For  the  cases  examined— upslope  propagation  over  both  a  lossless  and  lossy  bottom, 
both  URTEXT  and  BBCM  have  good  agreement  between  the  two.  In  the  close  range, 
BBCM  does  not  show  surface  interference  while  URTEXT  does.  However,  BBCM  does 
develop  a  pressure  field  in  the  underlying  bottom  while  URTEXT  is  effective  only  in  the 
water  column.  Furthermore,  BBCM  is  much  quicker  to  create  the  pressure  field.  Using  a 
range  resolution  of  25  m,  BBCM  created  a  pressure  field  for  each  scenario  in  2.5  hours  per 
run.  Using  the  same  equipment,  URTEXT  took  10  hours  each  run  to  develop  the  same 
resolution  without  calculating  transmission  loss  in  the  bottom.  BBCM  therefore  is  more 
time  effective. 

Other  strengths  of  BBCM  are  its  flexibility.  Because  the  input  files  model  of  the 
environment,  different  conditions  can  be  simulated  by  changing  the  input  files.  Upslope, 
downslope,  cross-canyon,  stratified  water  or  bottom  layers  can  all  be  simulated.  URTEXT  is 
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limited  only  to  the  condition  of  the  wedge  shaped  bottom.  It  is  still  a  good  tool  for  analysis 
of  the  specific  problem,  but  BBCM  is  a  more  flexible  tool  for  real  world  conditions. 

Areas  of  further  research  present  themselves.  URTEXT  needs  to  be  sped  up. 
MATLAB  is  very  fast  with  matrix  arithmetic  and  manipulation.  A  revision  of  URTEXT 
which  takes  advantage  of  this  capability  would  be  handy  in  the  analysis  of  the  wedge 
problem.  Also,  new  developments  to  the  image  problem  include  using  doublets  and  sets  of 
four  images.  Perhaps  including  these  as  subroutines  would  speed  up  URTEXT. 

The  examination  of  a  sloping  interface  over  a  rigid  sub-bottom  needs  to  be  investi gated. 
A  modification  of  URTEXT  to  include  this  condition  promises  interesting  results. 
Furthermore,  a  solution  of  the  above  problem  which  take  incident  angle  needs  to  be 
determined. 

Finally,  both  the  downslope  and  cross-slope  cases  need  to  be  investigated  with  BBCM. 
The  upslope  case  is  a  standard  comparison  with  other  models,  but  predictions  in  either 
conditions  might  prove  useful. 
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APPENDIX  A 

To  develop  a  comparison  standard  between  the  image  model  and  the  coupled  mode 
model,  we  translated  the  image  program  URTEXT  from  FORTRAN  to  MATLAB™  [Ref 
37].  MATLAB  is  a  high  level  scripting  language  which  runs  on  a  variety  of  computer 
operating  systems  including  VAX,  MS-DOS,  UNIX,  and  Macintosh.  This  improved  the 
flexibility  of  the  image  model  to  let  it  be  run  on  desktop  computers  and  advanced 
workstations.  Furthermore,  since  the  coupled  mode  model  is  also  written  in  MATLAB, 
comparison  of  CPU  run  times  can  be  made. 

Validation  of  the  MATLAB  version  of  URTEXT,  called  URTEXTBAT,  was  conducted 
by  comparing  the  results  of  URTEXTBAT  to  those  of  URTEXT  as  run  on  the  IBM  3033 
located  at  the  Naval  Postgraduate  School.  Table  A-l  shows  the  input  conditions  of  the 
programs.  Pressure  amplitudes  were  computed  for  scaling  distances  from  three  to  ten  dump 
distances  from  the  wedge  apex  and  at  receiver  angles  from  the  surface  to  the  bottom. 

TABLE  A-l 


p 

10° 

Y 

5° 

ri 

2 

P1/P2 

2.0 

ci/c2 

0.5 

a/k2 

0.001 

yo 

0.0 

Tables  A-2,  A-3,  and  A-4  show  the  results  of  the  comparison.  The  major  error  at  the 
surface  is  due  to  the  difference  in  precision  of  the  different  languages.  FORTRAN  uses 
single  precision  variables,  accurate  to  six  decimal  places.  MATLAB  is  precise  to  13  places, 
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therefore  at  the  small  pressures  near  the  surface  the  round-off  errors  are  extreme.  Excluding 
the  errors  at  the  boundaries,  the  average  absolute  difference  of  the  field  was  0.00206. 

TABLE  A-2 


FORTRAN                                         r2 

6 

3 

4 

5 

6 

7 

8 

9 

10 

10 

0.00324 

0.00813 

0.0026 

0.00158 

0.00118 

0.00118 

0.00153 

0.0031 

9 

1.1393 

0.75924 

0.63381 

0.4346 

0.30078 

0.2086 

0.13697 

0.0337 

8 

2.13245 

1.17088 

1.16923 

0.80956 

0.56091 

0.38328 

0.20159 

0.14456 

7 

2.86145 

1.4393 

1.52998 

1.07776 

0.74823 

0.47362 

0.34154 

0.34534 

6 
5 

3.24737 

1.6604 

1.66464 

1.20835 

0.84126 

0.5739 

0.56706 

0.5354 

3.26776 

1.93661 

1.56257 

1.19399 

0.8313 

0.79333 

0.74299 

0.67374 

4 
3 

2.97124 

2.06111 

1.25278 

1.05485 

0.91045 

0.9172 

0.82486 

0.72998 

2.51693 

1.97233 

0.80655 

0.84476 

0.98048 

0.90308 

0.79069 

0.68846 

2 
1 
0 

2.39419 

1.65642 

0.51872 

0.72407 

0.84396 

0.74538 

0.64007 

0.55031 

1.65191 

1.13991 

0.3643 

0.54176 

0.54487 

0.46663 

0.39299 

0.33296 

0.75925 

0.48278 

0.14244 

0.15337 

0.14283 

0.11185 

0.08627 

0.06728 

TABLE  A-3 


MATLAB™ 

r2 

6 

3 

4 

5 

6 

7 

8 

9 

10 

10 

9 

8 

0.0003 

0.0008 

0.0003 

0.0002 

0.0001 

0.0001 

0.0002 

0.0003 

1.1378 

0.7536 

0.6319 

0.4337 

0.3004 

0.2087 

0.1377 

0.0358 

2.1323 

1.1668 

1.1681 

0.8093 

0.5612 

0.3845 

0.2038 

0.1422 

7 

2.8627 

1.4368 

1.5293 

1.0781 

0.7494 

0.4822 

0.3388 

0.343 

6 

3.25 

1.6586 

1.6652 

1.2094 

0.8436 

0.571 

0.5643 

0.5336 

5 
4 

3.2716 

1.9359 

1.5637 

1.1955 

0.8286 

0.7897 

0.7412 

0.6726 

2.9753 

2.0614 

1.2543 

1.0564 

0.9016 

0.9152 

0.824 

0.7296 

3 

2.5161 

1.9736 

0.8081 

0.8451 

0.9778 

0.9026 

0.7908 

0.6887 

2 

2.3904 

1.6583 

0.5202 

0.721 

0.8439 

0.746 

0.6408 

0.551 

1 

1.6515 

1.1421 

0.3659 

0.5405 

0.5463 

0.468 

0.3941 

0.3339 

0 

0.7354 

0.485 

0.1428 

0.1557 

0.1449 

0.1134 

0.0875 

0.0683 

38 


TABLE  A-4 


Absolute  Difference           Ti 

6 

3 

4 

5 

6 

7 

8 

9 

10 

10 

0.00294 

0.00733 

0.00230 

0.00138 

0.00108 

0.00108 

0.00133 

0.00280 

9 
8 

0.001  50 

0.00564 

0.00191 

0.00090 

0.00038 

0.00010 

0.00073 

0.00210 

0.0001  5 

0.00408 

0.00113 

0.00026 

0.00029 

0.00122 

0.00221 

0.00236 

7 

0.00125 

0.00250 

0.00068 

0.00034 

0.00117 

0.00858 

0.00274 

0.00234 

6 

0.00263 

0.00180 

0.00056 

0.00105 

0.00234 

0.00290 

0.00276 

0.00180 

5 
4 

0.00384 

0.00071 

0.00113 

0.00151 

0.00270 

0.00363 

0.00179 

0.00114 

0.00406 

0.00029 

0.00152 

0.00155 

0.00885 

0.00200 

0.00086 

0.00038 

3 

0.00083 

0.00127 

0.00155 

0.00034 

0.00268 

0.00048 

0.0001 1 

0.00024 

2 

0.00379 

0.00188 

0.00148 

0.00307 

0.00006 

0.00062 

0.00073 

0.00069 

1 
0 

0.00041 

0.00219 

0.00160 

0.00126 

0.00143 

0.00137 

0.00111 

0.00094 

0.02385 

0.00222 

0.00036 

0.00233 

0.00207 

0.00155 

0.00123 

0.00102 

To  see  if  the  MATLAB  is  more  precise,  a  run  was  made  on  an  18°  slope  with  the  data 
taken  from  the  surface  to  2°  below  sampled  every  0.1°.  Figure  A-l  shows  pressure  falls  off 
linearly  as  the  sound  nears  the  surface.  This  implies  the  model  is  valid  at  the  boundaries. 
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Figure  A-l  -  Boundary  Condition  Check 

To  further  explore  the  differences  between  the  two  models,  a  third  model,  called 
WEDGES LO  was  used.  This  model  was  validated  by  Kaswandi  [Ref  14].  This  program  is 
written  in  BASIC  to  be  run  on  an  MS-DOS  based  computer,  but  used  double  precision 
variables.  A  single  run  at  r2=5  showed  correlation  between  WEDGES  LO  and 
URTEXTBAT  to  four  places,  including  at  the  boundary. 

Some  special  cases  were  run  on  URTEXTBAT  before  it  was  considered  valid.  First,  to 
verify  reciprocity  of  the  field,  the  source  was  set  at  ry=10  and  y  varied  from  0°  to  10°.  The 
receiver  was  fixed  at  6  =  5°  and  at  r^=2.  Table  A -5  shows  the  results  and  an  error 
comparison.  The  errors  are  within  the  round  off  tolerance  of  the  program,  showing 
URTEXT  as  valid  for  both  the  upslope  and  downslope  fields  as  well  as  validating  the 
acoustic  reciprocity  of  the  model. 
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TABLE  A-5  -  RECIPROCITY  CHECK 


6 

Pressure 

Difference 

10 

0 

0.00310 

9 

0.03349 

0.00021 

8 

0.14344 

0.00112 

7 

0.34385 

0.00149 

6 

0.53402 

0.00138 

5 

0.67274 

0.00100 

4 

0.72949 

0.00049 

3 

0.6885 

0.00004 

2 

0.55086 

0.00055 

1 

0.33387 

0.00091 

0 

0.06837 

0.00109 

Another  test  was  run  using  matched  acoustic  impedance  of  the  water  and  the  bottom. 
The  bottom  then  becomes  transparent  and  the  pressure  field  will  be  that  of  a  surface 
interference  effect  (Lloyd's  Mirror).  Density  and  speed  of  sound  ratios  were  set  to  .9999 
and  slope  angle  was  set  to  18°.  The  receiver  was  set  at  the  same  depth  as  the  source  and 
moved  downslope.  The  pressure  field  predicted  is  approximated  by  the  equation 

khd 


Hr.h)-  — 

r 


sin' 


(A-l) 


where  h  and  d  are  depths  of  the  receiver  and  the  source,  A  is  the  pressure  amplitude  at  the 
source  at  one  meter,  and  k  is  the  wave  number.  Frequency  used  was  1000  Hz. 

Figure  A -2  shows  a  plot  of  pressure  as  a  function  of  range  from  the  source. 
Interference  maxima  and  minima  from  URTEXTBAT  match  those  predicted  by  equation 
(A-l). 
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As  seen  from  the  above,  URTEXTBAT  has  undergone  extensive  proofing.  When 
compared  with  previous  versions  in  FORTRAN  and  BASIC,  it  is  more  flexible  both  in 
platforms  run  and  model  variables.  While  the  program  does  not  take  full  advantage  of  the 
features  of  MATLAB,  and  has  significant  room  for  optimizations,  the  program  in  its  present 
form  is  more  than  sufficient  for  the  purposes  of  this  research. 
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APPENDIX  B 

%%%%%%%%%%%%%%%%%%%%%%%% 
% 

%   URTEXT  for  MATLAB 

% 

%   Adapted  from  FORTRAN 

% 

%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Enter  Input  Variables 

%  Required  are:  Range  of  Source  and  Receiver 

%  Wedge  angle  in  Degrees 

%  Source  and  Reciever  angles  from  the  Bottom 

%  Speed  of  sound  ratio  and  density  ratios 

%  Bottom  loss  coefficient  and  cross  slope  range 

%%%%%%%%%%%%%%%%%%%%%%%% 

B=input( 'Enter  bottom  wedge  angle:  '); 
G=input( 'Enter  source  angle  from  the  surface:  '); 
D=input( 'Enter  receiver  angle  from  the  surface:  '); 
rhol=input( 'Enter  Density  ratio,  water  to  bottom:  '); 
CC=input( 'Enter  speed  ratio,  water  to  bottom:  '); 
rl=input( 'Enter  range  of  source  from  apex:  '); 
r2=input( 'Enter  range  of  receiver  from  apex:  '); 
AL=input( 'Enter  bottom  loss  coefficient:  '); 
yO=input( 'Enter  cross-slope  range:  '); 

%  Determine  number  of  image  pairs 

%  and  convert  to  radians.  Also  determine  constants 

%  Function  rads(A)  converts  degrees  into  radians 

%  Contained  as  a  seperate  M-File  in  the  MATLAB  Toolbox 

function  X=rads(A) 

%  Converts  degrees  into  radians 

X=A*pi/180; 

end 


Nl=fix(180/B) ; 
B=rads(B) ; 
G=rads(G) ; 
D=rads(D) ; 
c2=CC"2; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Determine  scaling  for  fast  or  slow  bottom 

%%%%%%%%%%%%%%%%%%%%%%%% 

tqq=tan(B) ; 
if  CC<1, 

tqql=acos( CC) ; 

tqq2=sin ( tqql ) ; 
else 

tqql=acos( 1/CC) ; 

tqq2=tan ( tqql ) ; 
end; 
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tqq3=2  * tqq2  *  tqq ; 

t4=pi/tqq3; 

ql=l/sqrt(2); 

d2=rl~2+r2~2+y0~2; 

r3=2*rl*r2; 

sl=l; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Determine  range  to  receiver  for  each  image 

%%%%%%%%%%%%%%%%%%%%%%%% 

for  n=l:l:Nl, 

if  sl>0 

tl(n)=(n-l)*B+G; 

else 

tl(n)=n*B-G; 

end; 

sl=-sl; 

r8(n)=sqrt(d2-r3*cos(tl(n)-D) ); 

r9(n)=sqrt(d2-r3*cos(tl(n)+D) )  ; 
end; 

pl=0;p2=0; 
for  n=l:l:Nl, 

s2=(-l)"(fix(n/2) ) ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calculate  Reflection  Coefficient  for  Upper  Space 
%%%%%%%%%%%%%%%%%%%%%%%% 

wl=2*c2*AL; 
il=fix( (n-l)/2); 
for  i=l:l:il, 

s(i)=abs(rl*sin(tl(n)-2*i*B)+r2*sin(2*i*B-D) )/r8(n) ; 

if  s(i)>=l, 
s(i)=l; 

end; 

c(i)=sqrt(l-(s(ir2)); 

t=s(i)/rhol; 

wO=(-c2+(c(i)A2)); 

y=sqrt( (wO*2)+(wlA2) ) ; 

z=abs(wO) ; 

if  y<=z, 
y=z; 

end; 

y l=ql *  sqrt ( y+wO ) ; 

y2=-ql*sqrt(y-w0) ; 

zl=t-y2; 

z2=-yl; 

z3=zl/(zl"2+z2~2); 

z4=-z2/(zl"2  +  z2's2); 

zl=t+y2; 

z2=yl; 

z5=zl*z3-z2*z4; 

z6=zl*z4+z2*z3; 

e(i)=z5; 

f(i)=z6; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Determine  Pressure  Contribution  for  each  upper  image 

%%%%%%%%%%%%%%%%%%%%%%%% 
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Zl=0;z2=0;z3=0;z4=0;z5=l;z6=0; 
if  n>2, 

for  i=l:l:il, 

zl=e(i) ; 

z2=f(i); 

z3=z5; 

z4=z6; 

z5=zl*z3-z2*z4; 

z6=zl*z4+z2*z3; 
end; 
end; 
zl=z5; 
z2=z6; 
t=t4*r8(n) ; 
z3=cos( t ) ; 
z4=-sin( t) ; 
z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
pl=pl+s2*z5/r8(n) ; 
p2=p2+s2*z6/r8(n) ; 
il=il+l; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Determine  Reflection  Coefficient  for  the  Lower  Images 
%%%%%%%%%%%%%%%%%%%%%%%% 

for  i=l:l:il, 

s(i)=abs(rl*sin(tl(n)-2*(i-l)*B)+r2*sin(2*(i-l)*B+D) )/r9(n) ; 
if  s(i)>l, 

s(i)=l; 
end; 

c(i)=sqrt(1.0001-s(i)~2); 
t=s(i) /rhol; 
w0=-c2+c(i) "2; 
y=sqrt(wO~2+wlA2) ; 
z=abs(wO) ; 
if  y<=z, 

y=z; 
end; 

yl=ql*sqrt(y+wO) ; 
y2=-ql *  sqrt ( y-wO ) ; 
zl=t-y2; 
z2=-yl ; 

z3=zl/(zl"2+z2'2); 
z4=-z2/(zl"2+z2"2); 
zl=t+y2; 
z2=yl; 

z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
e(i)=z5; 
f(i)=z6; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Find  Pressure  Contribution  for  Lower  Images 
%%%%%%%%%%%%%%%%%%%%%%%% 

zl=0 ; z2=0 ; z3=0 ; z4=0 ; z5=l ; z6=0 ; 
for  i=l:l:il, 
zl=e(i); 
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z2=f(i); 

z3=z5; 

z4=z6; 

z5=zl*z3-z2*z4 ; 

z6=zl*z4+z2*z3; 
end; 
zl=z5; 
z2=z6; 
t=t4*r9(n); 
z3=cos(t) ; 
z4=-sin( t) ; 
z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
pl=pl+s2*z5/r9(n) ; 
p2=p2+s2*z6/r9(n) ; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Total  Pressure  at  the  reciever  is: 

%%%%%%%%%%%%%%%%%%%%%%%% 

pz=sqrt(pl~2+p2~2) ; 
end; 
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%%%%%%%%%%%%%%%%%%%%%%%% 
% 

%     URTEXTBAT  for  MATLAB 

% 

%     Adapted  from  FORTRAN 

% 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  This  program  determines  a  2  dimensional 

%  Sound  pressure  field  over  a  sloping  bottom 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Enter  Input  Variables 

%  Required  are:  Range  of  Source 

%  Maximum  range  of  the  Receiver 

%  Wedge  angle  in  Degrees 

%  Source  angle  from  the  Surface 

%  Speed  of  sound  ratio  and  density  ratios 

%  Bottom  loss  coefficient  and  cross  slope  range. 

%  The  program  determines  receiver  angles  over  10  increments 

%  and  range  from  source  range  +  one  dump  distance  to  the 

%  maximum  receiver  range 

%%%%%%%%%%%%%%%%%%%%%%%% 

B=input( 'Enter  bottom  wedge  angle:  '  )  ; 

G=input( 'Enter  source  angle  from  the  surface:  '); 

rhol=input( 'Enter  Density  ratio,  water  to  bottom:  '); 

CC=input( 'Enter  speed  ratio,  water  to  bottom:  '); 

rl=input( 'Enter  range  of  source  from  apex:  '); 

range=input( 'Enter  maximum  receiver  range:  '); 

AL=input( 'Enter  bottom  loss  coefficient:  '); 

b=B/ 10 ; countk=0 ; count 1=0 ; 

y0=input( 'Enter  cross  slope  range:  '); 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Determine  number  of  image  pairs 

%  and  convert  angles  to  radians 

%%%%%%%%%%%%%%%%%%%%%%%% 

time=clock; 
Nl=fix(180/B); 
B=rads(B) ; 
G=rads(G) ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Loop  for  field  calculation 

%%%%%%%%%%%%%%%%%%%%%%%% 

for  k=rl+l :range/10 :range, 
r2=k; 

countk=countk+l ; 
for  l=0:b:b*10, 

D=l ; countl=countl+l ; 

D=rads(D) ; 

c2=CC~2; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Find  scaling  distance  for  fast  or  slow  bottom 

%%%%%%%%%%%%%%%%%%%%%%%% 

tqq=tan(B) ; 
if  CC<1, 
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tqql=acos(CC) ; 

tqq2=sin( tqql ) ; 
else 

tqql=acos ( 1 /CC ) ; 

tqq2=tan( tqql ) ; 
end; 

tqq3=2  * tqq2  *tqq; 
t4=pi/tqq3; 
ql=l/sqrt(2) ; 
d2=rl"2+r2~2+y0/v2; 
r3=2*rl*r2; 
■  1-1  j 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Determine  range  for  each  image  to  reciever 
%%%%%%%%%%%%%%%%%%%%%%%% 

for  n=l:l:Nl, 

if  sl>0 

tl(n)=(n-l)*B+G; 

else 

tl(n)=n*B-G; 

end; 

sl=-sl ; 

r8(n)=sqrt(d2-r3*cos(tl(n)-D) ) ; 

r9(n)=sqrt(d2-r3*cos(tl(n)+D) ) ; 
end; 

pl=0;p2=0; 
for  n=l:l:Nl, 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Find  Reflection  coefficient  for  upper  images 
%%%%%%%%%%%%%%%%%%%%%%%% 

s2=(-l)"(fix(n/2)); 
wl=2*c2*AL; 
il=fix((n-l)/2); 
for  i=l:l:il, 

s(i)=abs(rl*sin(tl(n)-2*i*B)+r2*sin(2*i*B-D) )/r8(n) ; 

if  s(i)>=l, 
s(i)=l; 

end; 

c(i)=sqrt(l-(s(ir2)); 

t=s( i ) /rhol ; 

w0=(-c2+(c(i)~2) ); 

y=sqrt( (w0"2 )+(wl"2 ) ) ; 

z=abs(w0 ) ; 

if  y<=z, 
y=z; 

end; 

yl=ql*sqrt(y+wO) ; 

y2=-ql*sqrt(y-w0) ; 

zl=t-y2; 

z2=-yl; 

z3=zl/(zl'v2  +  z2~2); 

z4=-z2/(zl"2+z2/v2); 

zl=t+y2; 

z2=yl; 

z5=zl*z3-z2*z4; 

z6=zl*z4+z2*z3; 

e(i)=z5; 
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f(i)=z6; 


f(i)=z6; 
end; 

zl=0; z2=0; z3=0; z4=0; z5=l; z6=0; 
if  n>2, 

for  i=l:l:il, 


-z2*z4; 
kz4+z2*z3; 


%%%%%%%%%%%%%%%%%%%%%%%% 

%  Find  Pressure  contribution  of  upper  images 

%%%%%%%%%%%%%%%%%%%%%%%% 

zl=z5; 

z2=z6; 

t=t4*r8(n) ; 

z3=cos( t) ; 

z4=-sin( t) ; 

z5=zl*z3-z2*z4; 

z6=zl*z4+z2*z3; 

pl=pl+s2*z5/r8(n); 

p2=p2+s2*z6/r8(n) ; 

il=il+l; 

for  i=l:l:il, 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Find  refection  coefficient  for  lower  images 
%%%%%%%%%%%%%%%%%%%%%%%% 

s(i)=abs(rl*sin(tl(n)-2*(i-l)*B)+r2*sin(2*(i-l)*B+D) )/r9(n) ; 

if  s(i)>l, 
s(i)=l; 

end; 

c(i)=sqrt(1.0001-s(ip2)  ; 

t=s( i ) /rhol ; 

w0=-c2+c(i) "2; 

y=sqrt(w0~2+wl"2) ; 

z=abs( wO ) ; 

if  y<=z, 
y=z; 

end; 

yl=ql*sqrt (y+wO) ; 

y2=-ql*sqrt(y-w0) ; 

zl=t-y2; 

z2=-yl; 

z3=zl/(zl~2+z2~2) ; 

z4=-z2/(zl~2+z2~2); 

zl=t+y2; 

z2=yl; 

z5=zl*z3-z2*z4; 

z6=zl*z4+z2*z3; 

e(i)=z5; 

f(i)=z6; 
end; 

zl=0; z2=0;z3=0;z4=0;z5=l; z6=0; 
for  i=l:l:il, 
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zl=e(i); 
z2=f(i); 
z3=z5; 
z4=z6; 

z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Find  Pressure  contribution  for  lower  images 
%%%%%%%%%%%%%%%%%%%%%%%% 

zl=z5; 
z2=z6; 
t=t4*r9(n); 
z3=cos( t ) ; 
z4=-sin( t ) ; 
z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
pl=pl+s2*z5/r9(n); 
p2=p2+s2*z6/r9(n) ; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Total  pressure  at  the  reciever  is: 
%%%%%%%%%%%%%%%%%%%%%%%% 

pz(countl,countk)=sqrt(pl~2+p2~2 ) ; 

%%%%%%%%%%%%%%%%%%%%%%%% 
%  Reset  variables 
%%%%%%%%%%%%%%%%%%%%%%%% 

clear  c; clear  c2; clear  d2; clear  e;  clear  f; clear  i;  clear  il; 

clear  n;clear  pl;clear  p2;clear  ql;clear  r3;clear  r8;  clear  r9; 

clear  s;clear  sl;clear  s2;clear  t;clear  tl; 

clear  tqq; clear  tqql; clear  tqq2; clear  tqq3; clear  wO;  clear  wl ; 

clear  y;clear  yl;clear  y2;clear  z;clear  zl;clear  z2;  clear  z3; 

clear  z4; clear  z5; clear  z6; 

end; 
countl=0; 
et=etime (clock , time ) 
end; 
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%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  SCENEGEN.M 

% 

%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  This  program  creates  input  files  for  BBCM, 

%  the  coupled  mode  model  by  Dr.  C.S.  Chiu. 

%  This  program  creates  the  inputs  for  a  wedge  shaped 

%  ocean  overlying  a  fluid  bottom. 

%  Both  water  and  bottom  are  considered  uniform  in  characteristics 

%  in  depth.  No  irregularities  in  the  bottom  depth  are 

%  generated,  nor  stratifications  in  sound  velocity  profile. 

% 

%%%%%%%%%%%%%%%%%%%%%%%% 

%  Input  all  the  pertinent  data 

B=input( "Enter  wedge  angle:  '); 

B=rads(B) ; 

z=input( 'Enter  maximum  depth  of  the  water  column:  '); 

w=input( 'Enter  length  of  the  waveguide:  '); 

dw=input( 'Enter  range  resolution:  '); 

dz=input( 'Enter  depth  resolution:  '); 

cl=input( 'What  is  the  speed  of  sound  in  the  water?  ' ) ; 

c2=input( 'What  is  the  speed  of  sound  in  the  bottom?  '); 

rl=input( 'What  is  the  density  of  the  water?  '); 

r2=input( 'What  is  the  density  of  the  bottom?  '); 

n=input( 'How  many  modes?  ' ) ; 

ak=input( 'What  is  the  bottom  loss  coefficient?'); 

f=input( 'What  is  the  source  frequency?  '); 

%  Build  data  files  —  range 

r=0 :dw:w; 

save  r.dat  r  /ascii; 

%  receiver  depth 

zr=0:10:z; 

save  zr.dat  zr  /ascii; 

%  attenuation  matrix 
for  i=l:l:n, 

attenu(n)=ak; 
end; 
save  attenu.dat  attenu  /ascii; 

%  sound  speed  profile 
a=w/dw;  b=z/dz; 
t=tan(B) ; 
for  i=l:l:a+l, 
for  j=l:l:b+l, 

if  t>( j-l)*dz/(w-(i-l)*dw) , 

c( j,i)=cl; 
else 

c ( j , i ) =c2 ; 
end; 
end; 
end; 

dummy =c( 1 :b, : ) ; 
for  i=l:l:b, 

dummy ( i , 1 ) =c2 ; 
end; 
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for  i=2:l:a, 

dummy  (  : ,  i )  =dummy  (  : ,  i- 1 ) ; 
end; 

c=  [  c ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy  ] ; 
save  c.dat  c  /ascii; 

%  density  profile 
for  i=l:l:a+l, 
for  j=l:l:b+l, 

if  t>( j-l)*dz/(w-(i-l)*dw), 

rho( j,i)=rl; 
else 

rho( j,i)=r2; 
end; 
end; 
end; 

dummy=rho( 1 :b, : ) ; 
for  i=l:l:b, 

dummy (i, 1 )=r2; 
end; 
for  i=2:l:a, 

dummy ( : ,i)=dummy( : ,i-l ) ; 
end; 

r ho=  [  r ho ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy ;  dummy  ] ; 
save  rho.dat  rho  /ascii; 

%  Source  Spectrum  Density  Profile 

ssp=[0  0;f-0.005  l;f  l;f+0.005  1;1000  0 ] ; 

save  ssp.dat  ssp  /ascii; 
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